Survival Analysis

Colorectal Cancer | Pyroptosis | Unique Genes per RCD Type | Gene Expression of Normal Samples

I. Preliminaries

Loading libraries

library("tidyverse")
library("tibble")
library("msigdbr")
library("ggplot2")
library("TCGAbiolinks")
library("RNAseqQC")
library("DESeq2")
library("ensembldb")
library("purrr")
library("magrittr")
library("vsn")
library("matrixStats")
library("dplyr")
library("grex")
library("survminer")
library("survival")

II. Downloading the TCGA gene expression data

Create a function for downloading TCGA gene expression data.

For more detailed documentation, refer to 2. Differential Gene Expression Analysis - TCGA.Rmd.

GDC_DIR = "../data/public/GDCdata"

query_and_filter_samples <- function(project) {
  query_tumor <- GDCquery(
    project = project,
    data.category = "Transcriptome Profiling",
    data.type = "Gene Expression Quantification",
    experimental.strategy = "RNA-Seq",
    workflow.type = "STAR - Counts",
    access = "open",
    sample.type = "Primary Tumor"
  )
  tumor <- getResults(query_tumor)

  query_normal <- GDCquery(
    project = project,
    data.category = "Transcriptome Profiling",
    data.type = "Gene Expression Quantification",
    experimental.strategy = "RNA-Seq",
    workflow.type = "STAR - Counts",
    access = "open",
    sample.type = "Solid Tissue Normal"
  )
  normal <- getResults(query_normal)

  submitter_ids <- inner_join(tumor, normal, by = "cases.submitter_id") %>%
    dplyr::select(cases.submitter_id)
  tumor <- tumor %>%
    dplyr::filter(cases.submitter_id %in% submitter_ids$cases.submitter_id)
  normal <- normal %>%
    dplyr::filter(cases.submitter_id %in% submitter_ids$cases.submitter_id)

  samples <- rbind(tumor, normal)
  unique(samples$sample_type)

  query_project <- GDCquery(
    project = project,
    data.category = "Transcriptome Profiling",
    data.type = "Gene Expression Quantification",
    experimental.strategy = "RNA-Seq",
    workflow.type = "STAR - Counts",
    access = "open",
    sample.type = c("Solid Tissue Normal", "Primary Tumor"),
    barcode = as.list(samples$sample.submitter_id)
  )

  # If this is your first time running this notebook (i.e., you have not yet downloaded the results of the query in the previous block),
  # uncomment the code block below

  # GDCdownload(
  #   query_coad,
  #   directory = GDC_DIR
  # )

  return(list(samples = samples, query_project = query_project))
}

Download the TCGA gene expression data for colorectal cancer (TCGA-COAD).

projects <- c("TCGA-COAD")

with_results_projects <- c()

samples <- list()
project_data <- list()

for (project in projects) {
  result <- tryCatch(
    {
      result <- query_and_filter_samples(project)
      samples[[project]] <- result$samples
      project_data[[project]] <- result$query_project

      with_results_projects <- c(with_results_projects, project)
    },
    error = function(e) {

    }
  )
}

Running the code block above should generate and populate a directory named GDCdata.

III. Data preprocessing

Construct the RNA-seq count matrix for each cancer type.

tcga_data <- list()
tcga_matrix <- list()

projects <- with_results_projects
for (project in projects) {
  tcga_data[[project]] <- GDCprepare(
    project_data[[project]], 
    directory = GDC_DIR,
    summarizedExperiment = TRUE
  )
}

for (project in projects) {
  count_matrix <- assay(tcga_data[[project]], "unstranded")

  # Remove duplicate entries
  count_matrix_df <- data.frame(count_matrix)
  count_matrix_df <- count_matrix_df[!duplicated(count_matrix_df), ]
  count_matrix <- data.matrix(count_matrix_df)
  rownames(count_matrix) <- cleanid(rownames(count_matrix))
  count_matrix <- count_matrix[!(duplicated(rownames(count_matrix)) | duplicated(rownames(count_matrix), fromLast = TRUE)), ]

  tcga_matrix[[project]] <- count_matrix
}

Format the samples table so that it can be fed as input to DESeq2.

for (project in projects) {
  rownames(samples[[project]]) <- samples[[project]]$cases
  samples[[project]] <- samples[[project]] %>%
    dplyr::select(case = "cases.submitter_id", type = "sample_type")
  samples[[project]]$type <- str_replace(samples[[project]]$type, "Solid Tissue Normal", "normal")
  samples[[project]]$type <- str_replace(samples[[project]]$type, "Primary Tumor", "tumor")
}

DESeq2 requires the row names of samples should be identical to the column names of count_matrix.

for (project in projects) {
  colnames(tcga_matrix[[project]]) <- gsub(x = colnames(tcga_matrix[[project]]), pattern = "\\.", replacement = "-")
  tcga_matrix[[project]] <- tcga_matrix[[project]][, rownames(samples[[project]])]

  # Sanity check
  print(all(colnames(tcga_matrix[[project]]) == rownames(samples[[project]])))
}

IV. Differential gene expression analysis

For more detailed documentation on obtaining the gene set, refer to 7. Differential Gene Expression Analysis - TCGA - Pan-cancer - Unique Genes.Rmd.

RCDdb <- "../data/public/rcd-gene-list/unique-genes/necroptosis-ferroptosis-pyroptosis/"

Write utility functions for filtering the gene sets, performing differential gene expression analysis, plotting the results, and performing variance-stabilizing transformation.

filter_gene_set_and_perform_dgea <- function(genes) {
  tcga_rcd <- list()

  for (project in projects) {
    rownames(genes) <- genes$gene_id
    tcga_rcd[[project]] <- tcga_matrix[[project]][rownames(tcga_matrix[[project]]) %in% genes$gene_id, ]
    tcga_rcd[[project]] <- tcga_rcd[[project]][, rownames(samples[[project]])]
  }

  dds_rcd <- list()
  res_rcd <- list()

  for (project in projects) {
    print(project)
    print("=============")
    dds <- DESeqDataSetFromMatrix(
      countData = tcga_rcd[[project]],
      colData = samples[[project]],
      design = ~type
    )
    dds <- filter_genes(dds, min_count = 10)
    dds$type <- relevel(dds$type, ref = "normal")
    dds_rcd[[project]] <- DESeq(dds)
    res_rcd[[project]] <- results(dds_rcd[[project]])
  }

  deseq.bbl.data <- list()

  for (project in projects) {
    deseq.results <- res_rcd[[project]]
    deseq.bbl.data[[project]] <- data.frame(
      row.names = rownames(deseq.results),
      baseMean = deseq.results$baseMean,
      log2FoldChange = deseq.results$log2FoldChange,
      lfcSE = deseq.results$lfcSE,
      stat = deseq.results$stat,
      pvalue = deseq.results$pvalue,
      padj = deseq.results$padj,
      cancer_type = project,
      gene_symbol = genes[rownames(deseq.results), "gene"]
    )
  }

  deseq.bbl.data.combined <- bind_rows(deseq.bbl.data)
  deseq.bbl.data.combined <- dplyr::filter(deseq.bbl.data.combined, abs(log2FoldChange) >= 1.5 & padj < 0.05)

  return(deseq.bbl.data.combined)
}

plot_dgea <- function(deseq.bbl.data.combined) {
  sizes <- c("<10^-15" = 4, "10^-10" = 3, "10^-5" = 2, "0.05" = 1)

  deseq.bbl.data.combined <- deseq.bbl.data.combined %>%
    mutate(fdr_category = cut(padj,
      breaks = c(-Inf, 1e-15, 1e-10, 1e-5, 0.05),
      labels = c("<10^-15", "10^-10", "10^-5", "0.05"),
      right = FALSE
    ))

  top_genes <- deseq.bbl.data.combined %>%
    group_by(cancer_type) %>%
    mutate(rank = rank(-abs(log2FoldChange))) %>%
    dplyr::filter(rank <= 10) %>%
    ungroup()

  ggplot(top_genes, aes(y = cancer_type, x = gene_symbol, size = fdr_category, fill = log2FoldChange)) +
    geom_point(alpha = 0.5, shape = 21, color = "black") +
    scale_size_manual(values = sizes) +
    scale_fill_gradient2(low = "blue", mid = "white", high = "red", limits = c(min(deseq.bbl.data.combined$log2FoldChange), max(deseq.bbl.data.combined$log2FoldChange))) +
    theme_minimal() +
    theme(
      axis.text.x = element_text(size = 9, angle = 90, hjust = 1)
    ) +
    theme(legend.position = "bottom") +
    theme(legend.position = "bottom") +
    labs(size = "Adjusted p-value", fill = "log2 FC", y = "Cancer type", x = "Gene")
}

perform_vsd <- function(genes) {
  tcga_rcd <- list()

  for (project in projects) {
    rownames(genes) <- genes$gene_id
    tcga_rcd[[project]] <- tcga_matrix[[project]][rownames(tcga_matrix[[project]]) %in% genes$gene_id, ]
    tcga_rcd[[project]] <- tcga_rcd[[project]][, rownames(samples[[project]])]
  }

  vsd_rcd <- list()

  for (project in projects) {
    print(project)
    print("=============")
    dds <- DESeqDataSetFromMatrix(
      countData = tcga_rcd[[project]],
      colData = samples[[project]],
      design = ~type
    )
    dds <- filter_genes(dds, min_count = 10)

    # Perform variance stabilization
    dds <- estimateSizeFactors(dds)
    nsub <- sum(rowMeans(counts(dds, normalized = TRUE)) > 10)
    vsd <- vst(dds, nsub = nsub)
    vsd_rcd[[project]] <- assay(vsd)
  }

  return(vsd_rcd)
}

Pyroptosis

Fetch the gene set of interest.

genes <- read.csv(paste0(RCDdb, "Pyroptosis.csv"))
print(genes)

genes$gene_id <- cleanid(genes$gene_id)
genes <- distinct(genes, gene_id, .keep_all = TRUE)
genes <- subset(genes, gene_id != "")
genes

Filter the genes to include only those in the gene set of interest, and then perform differential gene expression analysis.

deseq.bbl.data.combined <- filter_gene_set_and_perform_dgea(genes)

[1] "TCGA-COAD"
[1] "============="

Warning: some variables in design formula are characters, converting to factorsestimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
-- replacing outliers and refitting for 3 genes
-- DESeq argument 'minReplicatesForReplace' = 7 
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing

deseq.bbl.data.combined

Plot the results.

plot_dgea(deseq.bbl.data.combined)

Perform variance-stabilizing transformation for further downstream analysis (i.e., for survival analysis).

vsd <- perform_vsd(genes)

[1] "TCGA-COAD"
[1] "============="

V. Downloading the clinical data

Download clinical data from TCGA, and perform some preprocessing: - The deceased column should be FALSE if the patient is alive and TRUE otherwise - The overall_survival column should reflect the follow-up time if the patient is alive and the days to death otherwise

download_clinical_data <- function(project) {
  clinical_data <- GDCquery_clinic(project)
  clinical_data$deceased <- ifelse(clinical_data$vital_status == "Alive", FALSE, TRUE)
  clinical_data$overall_survival <- ifelse(clinical_data$vital_status == "Alive",
    clinical_data$days_to_last_follow_up,
    clinical_data$days_to_death
  )

  return(clinical_data)
}

tcga_clinical <- list()
for (project in projects) {
  tcga_clinical[[project]] <- download_clinical_data(project)
}

VI. Performing survival analysis

Write utility functions for performing survival analysis.

construct_gene_df <- function(gene_of_interest, project) {
  gene_df <- vsd[[project]] %>%
    as.data.frame() %>%
    rownames_to_column(var = "gene_id") %>%
    gather(key = "case_id", value = "counts", -gene_id) %>%
    left_join(., genes, by = "gene_id") %>%
    dplyr::filter(gene == gene_of_interest) %>%
    dplyr::filter(case_id %in% rownames(samples[[project]] %>% dplyr::filter(type == "normal")))

  q1 <- quantile(gene_df$counts, probs = 0.25)
  q3 <- quantile(gene_df$counts, probs = 0.75)
  gene_df$strata <- ifelse(gene_df$counts >= q3, "HIGH", ifelse(gene_df$counts <= q1, "LOW", "MIDDLE"))
  gene_df <- gene_df %>% dplyr::filter(strata %in% c("LOW", "HIGH"))
  gene_df$case_id <- paste0(sapply(strsplit(as.character(gene_df$case_id), "-"), `[`, 1), '-',
                          sapply(strsplit(as.character(gene_df$case_id), "-"), `[`, 2), '-', 
                          sapply(strsplit(as.character(gene_df$case_id), "-"), `[`, 3))
  gene_df <- merge(gene_df, tcga_clinical[[project]], by.x = "case_id", by.y = "submitter_id")
  
  return(gene_df)
}

compute_surival_fit <- function(gene_df) {
  return (survfit(Surv(overall_survival, deceased) ~ strata, data = gene_df))
}

compute_cox <- function(gene_df) {
  return (coxph(Surv(overall_survival, deceased) ~ strata, data=gene_df))
}

plot_survival <- function(fit) {
  return(ggsurvplot(fit,
    data = gene_df,
    pval = T,
    risk.table = T,
    risk.table.height = 0.3
  ))
}

compute_survival_diff <- function(gene_df) {
  return(survdiff(Surv(overall_survival, deceased) ~ strata, data = gene_df))
}

Perform survival analysis by testing for the difference in the Kaplan-Meier curves using the G-rho family of Harrington and Fleming tests: https://rdrr.io/cran/survival/man/survdiff.html

Our genes of interest are GSDMD (the primary executor of pyroptosis) and the differentially expressed genes.

significant_projects <- c()
significant_genes <- c()

ctr <- 1
for (project in projects) {
  for (gene in c("GSDMD", genes$gene)) {
    cat(project, gene, "\n\n")
    error <- tryCatch (
      {
        gene_df <- construct_gene_df(gene, project)
      },
      error = function(e) {
        cat("\n\n============================\n\n")
        e
      }
    )
    
    if(inherits(error, "error")) next

    if (nrow(gene_df) > 0) {
      fit <- compute_surival_fit(gene_df)
      tryCatch (
        {
          survival <- compute_survival_diff(gene_df)
          cox <- compute_cox(gene_df)
          print(ctr)
          ctr <- ctr + 1
          print(survival)
          cat("\n")
          print(cox)
          print(plot_survival(fit))
          if (pchisq(survival$chisq, length(survival$n)-1, lower.tail = FALSE) < 0.05) {
            significant_projects <- c(significant_projects, project)
            significant_genes <- c(significant_genes, gene)
          }
        },
        error = function(e) {
        }
      )
      
    }
    
    cat("\n\n============================\n\n")
  }
}

TCGA-COAD GSDMD 

[1] 1
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=5, 17 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 3        3     1.93     0.589      1.18
strata=LOW  2        2     3.07     0.371      1.18

 Chisq= 1.2  on 1 degrees of freedom, p= 0.3 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)     z     p
strataLOW -1.2084    0.2987   1.1730 -1.03 0.303

Likelihood ratio test=1.23  on 1 df, p=0.2675
n= 5, number of events= 5 
   (17 observations deleted due to missingness)


============================

TCGA-COAD CHMP7 

[1] 2
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=7, 20 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 4        4     4.59    0.0754     0.248
strata=LOW  3        3     2.41    0.1434     0.248

 Chisq= 0.2  on 1 degrees of freedom, p= 0.6 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

            coef exp(coef) se(coef)     z     p
strataLOW 0.4115    1.5090   0.8315 0.495 0.621

Likelihood ratio test=0.24  on 1 df, p=0.6217
n= 7, number of events= 7 
   (20 observations deleted due to missingness)


============================

TCGA-COAD GSDMC 

[1] 3
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=6, 16 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 4        4     4.88     0.160     0.985
strata=LOW  2        2     1.12     0.699     0.985

 Chisq= 1  on 1 degrees of freedom, p= 0.3 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

            coef exp(coef) se(coef)     z     p
strataLOW 0.9671    2.6303   1.0107 0.957 0.339

Likelihood ratio test=0.88  on 1 df, p=0.3471
n= 6, number of events= 6 
   (16 observations deleted due to missingness)


============================

TCGA-COAD ELANE 

[1] 4
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=8, 14 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 5        5     3.76     0.407     0.904
strata=LOW  3        3     4.24     0.361     0.904

 Chisq= 0.9  on 1 degrees of freedom, p= 0.3 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)      z     p
strataLOW -0.7860    0.4557   0.8464 -0.929 0.353

Likelihood ratio test=0.94  on 1 df, p=0.3319
n= 8, number of events= 8 
   (14 observations deleted due to missingness)


============================

TCGA-COAD IRF1 

[1] 5
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=5, 17 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 3        3     4.35     0.419      4.26
strata=LOW  2        2     0.65     2.804      4.26

 Chisq= 4.3  on 1 degrees of freedom, p= 0.04 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

               coef exp(coef)  se(coef)     z     p
strataLOW 2.232e+01 4.923e+09 3.308e+04 0.001 0.999

Likelihood ratio test=4.61  on 1 df, p=0.03188
n= 5, number of events= 5 
   (17 observations deleted due to missingness)


============================

TCGA-COAD CYCS 

[1] 6
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=6, 16 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 1        1     1.45    0.1397     0.211
strata=LOW  5        5     4.55    0.0445     0.211

 Chisq= 0.2  on 1 degrees of freedom, p= 0.6 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

            coef exp(coef) se(coef)     z     p
strataLOW 0.5241    1.6889   1.1513 0.455 0.649

Likelihood ratio test=0.23  on 1 df, p=0.6342
n= 6, number of events= 6 
   (16 observations deleted due to missingness)


============================

TCGA-COAD GSDMA 

[1] 7
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=8, 14 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 3        3     2.54    0.0845     0.149
strata=LOW  5        5     5.46    0.0393     0.149

 Chisq= 0.1  on 1 degrees of freedom, p= 0.7 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)      z     p
strataLOW -0.3162    0.7289   0.8231 -0.384 0.701

Likelihood ratio test=0.15  on 1 df, p=0.7014
n= 8, number of events= 8 
   (14 observations deleted due to missingness)


============================

TCGA-COAD CASP4 

[1] 8
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=7, 15 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 2        2    0.902     1.335      1.74
strata=LOW  5        5    6.098     0.198      1.74

 Chisq= 1.7  on 1 degrees of freedom, p= 0.2 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)      z     p
strataLOW -1.2525    0.2858   1.0088 -1.242 0.214

Likelihood ratio test=1.45  on 1 df, p=0.2278
n= 7, number of events= 7 
   (15 observations deleted due to missingness)


============================

TCGA-COAD BAK1 

[1] 9
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=7, 15 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 2        2      2.1   0.00499    0.0079
strata=LOW  5        5      4.9   0.00214    0.0079

 Chisq= 0  on 1 degrees of freedom, p= 0.9 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

            coef exp(coef) se(coef)     z     p
strataLOW 0.0781    1.0812   0.8787 0.089 0.929

Likelihood ratio test=0.01  on 1 df, p=0.9289
n= 7, number of events= 7 
   (15 observations deleted due to missingness)


============================

TCGA-COAD NOD1 

[1] 10
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=6, 16 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 4        4      3.6    0.0444     0.135
strata=LOW  2        2      2.4    0.0667     0.135

 Chisq= 0.1  on 1 degrees of freedom, p= 0.7 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)      z     p
strataLOW -0.3429    0.7097   0.9368 -0.366 0.714

Likelihood ratio test=0.14  on 1 df, p=0.712
n= 6, number of events= 6 
   (16 observations deleted due to missingness)


============================

TCGA-COAD NLRP7 

[1] 11
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=6, 16 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 2        2    0.983     1.051      1.59
strata=LOW  4        4    5.017     0.206      1.59

 Chisq= 1.6  on 1 degrees of freedom, p= 0.2 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)      z     p
strataLOW -1.4388    0.2372   1.2359 -1.164 0.244

Likelihood ratio test=1.46  on 1 df, p=0.2262
n= 6, number of events= 6 
   (16 observations deleted due to missingness)


============================

TCGA-COAD CASP3 

[1] 12
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=7, 15 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 2        2      1.6    0.0987     0.147
strata=LOW  5        5      5.4    0.0293     0.147

 Chisq= 0.1  on 1 degrees of freedom, p= 0.7 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)      z     p
strataLOW -0.3505    0.7044   0.9188 -0.381 0.703

Likelihood ratio test=0.14  on 1 df, p=0.7067
n= 7, number of events= 7 
   (15 observations deleted due to missingness)


============================

TCGA-COAD GSDMB 

[1] 13
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=3, 19 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 2        2    2.667     0.167         2
strata=LOW  1        1    0.333     1.333         2

 Chisq= 2  on 1 degrees of freedom, p= 0.2 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

               coef exp(coef)  se(coef) z p
strataLOW 2.215e+01 4.171e+09 4.566e+04 0 1

Likelihood ratio test=2.2  on 1 df, p=0.1383
n= 3, number of events= 3 
   (19 observations deleted due to missingness)


============================

TCGA-COAD GZMB 

[1] 14
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=7, 15 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 3        3     1.91     0.619      1.03
strata=LOW  4        4     5.09     0.233      1.03

 Chisq= 1  on 1 degrees of freedom, p= 0.3 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)      z     p
strataLOW -0.9149    0.4005   0.9283 -0.986 0.324

Likelihood ratio test=1  on 1 df, p=0.3178
n= 7, number of events= 7 
   (15 observations deleted due to missingness)


============================

TCGA-COAD GSDME 

[1] 15
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=6, 16 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 4        4    5.017     0.206      1.59
strata=LOW  2        2    0.983     1.051      1.59

 Chisq= 1.6  on 1 degrees of freedom, p= 0.2 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

           coef exp(coef) se(coef)     z     p
strataLOW 1.439     4.215    1.236 1.164 0.244

Likelihood ratio test=1.46  on 1 df, p=0.2262
n= 6, number of events= 6 
   (16 observations deleted due to missingness)


============================

TCGA-COAD CHMP3 

[1] 16
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=7, 15 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 3        3     2.25     0.254     0.417
strata=LOW  4        4     4.75     0.120     0.417

 Chisq= 0.4  on 1 degrees of freedom, p= 0.5 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)      z     p
strataLOW -0.5293    0.5890   0.8287 -0.639 0.523

Likelihood ratio test=0.4  on 1 df, p=0.5252
n= 7, number of events= 7 
   (15 observations deleted due to missingness)


============================

TCGA-COAD DPP9 

[1] 17
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=8, 14 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 6        6     6.35    0.0190     0.102
strata=LOW  2        2     1.65    0.0731     0.102

 Chisq= 0.1  on 1 degrees of freedom, p= 0.7 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

            coef exp(coef) se(coef)     z    p
strataLOW 0.2777    1.3200   0.8717 0.319 0.75

Likelihood ratio test=0.1  on 1 df, p=0.754
n= 8, number of events= 8 
   (14 observations deleted due to missingness)


============================

TCGA-COAD NOD2 

[1] 18
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=6, 16 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 3        3     1.73     0.926      1.67
strata=LOW  3        3     4.27     0.376      1.67

 Chisq= 1.7  on 1 degrees of freedom, p= 0.2 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

            coef exp(coef) se(coef)      z     p
strataLOW -1.395     0.248    1.164 -1.198 0.231

Likelihood ratio test=1.69  on 1 df, p=0.1937
n= 6, number of events= 6 
   (16 observations deleted due to missingness)


============================

TCGA-COAD NLRC4 

[1] 19
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=7, 15 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 3        3     1.41     1.786      2.82
strata=LOW  4        4     5.59     0.451      2.82

 Chisq= 2.8  on 1 degrees of freedom, p= 0.09 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)      z     p
strataLOW -1.7438    0.1748   1.1660 -1.496 0.135

Likelihood ratio test=2.68  on 1 df, p=0.1015
n= 7, number of events= 7 
   (15 observations deleted due to missingness)


============================

TCGA-COAD GSDMD 

[1] 20
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=5, 17 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 3        3     1.93     0.589      1.18
strata=LOW  2        2     3.07     0.371      1.18

 Chisq= 1.2  on 1 degrees of freedom, p= 0.3 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)     z     p
strataLOW -1.2084    0.2987   1.1730 -1.03 0.303

Likelihood ratio test=1.23  on 1 df, p=0.2675
n= 5, number of events= 5 
   (17 observations deleted due to missingness)


============================

TCGA-COAD TIRAP 

[1] 21
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=5, 17 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 2        2     0.65     2.804      4.26
strata=LOW  3        3     4.35     0.419      4.26

 Chisq= 4.3  on 1 degrees of freedom, p= 0.04 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

                coef  exp(coef)   se(coef)      z     p
strataLOW -2.232e+01  2.031e-10  3.308e+04 -0.001 0.999

Likelihood ratio test=4.61  on 1 df, p=0.03188
n= 5, number of events= 5 
   (17 observations deleted due to missingness)


============================

TCGA-COAD SCAF11 

[1] 22
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=7, 15 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 4        4      3.8    0.0100    0.0262
strata=LOW  3        3      3.2    0.0119    0.0262

 Chisq= 0  on 1 degrees of freedom, p= 0.9 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)      z     p
strataLOW -0.1336    0.8749   0.8268 -0.162 0.872

Likelihood ratio test=0.03  on 1 df, p=0.8716
n= 7, number of events= 7 
   (15 observations deleted due to missingness)


============================

TCGA-COAD NLRP6 

[1] 23
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=8, 14 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 3        3     2.87   0.00587    0.0101
strata=LOW  5        5     5.13   0.00328    0.0101

 Chisq= 0  on 1 degrees of freedom, p= 0.9 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

              coef exp(coef) se(coef)      z    p
strataLOW -0.07751   0.92542  0.77060 -0.101 0.92

Likelihood ratio test=0.01  on 1 df, p=0.92
n= 8, number of events= 8 
   (14 observations deleted due to missingness)


============================

TCGA-COAD AIM2 

[1] 24
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=8, 14 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 3        3     4.04     0.266     0.636
strata=LOW  5        5     3.96     0.271     0.636

 Chisq= 0.6  on 1 degrees of freedom, p= 0.4 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

            coef exp(coef) se(coef)     z     p
strataLOW 0.6599    1.9347   0.8422 0.784 0.433

Likelihood ratio test=0.66  on 1 df, p=0.4153
n= 8, number of events= 8 
   (14 observations deleted due to missingness)


============================

TCGA-COAD CASP6 

[1] 25
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=8, 14 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 2        2     1.85   0.01176    0.0172
strata=LOW  6        6     6.15   0.00354    0.0172

 Chisq= 0  on 1 degrees of freedom, p= 0.9 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)      z     p
strataLOW -0.1145    0.8918   0.8720 -0.131 0.896

Likelihood ratio test=0.02  on 1 df, p=0.8962
n= 8, number of events= 8 
   (14 observations deleted due to missingness)


============================

TCGA-COAD NLRP2 

[1] 26
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=7, 15 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 2        2     3.35     0.546      1.32
strata=LOW  5        5     3.65     0.501      1.32

 Chisq= 1.3  on 1 degrees of freedom, p= 0.3 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

           coef exp(coef) se(coef)     z     p
strataLOW 1.199     3.318    1.105 1.085 0.278

Likelihood ratio test=1.48  on 1 df, p=0.2236
n= 7, number of events= 7 
   (15 observations deleted due to missingness)


============================

TCGA-COAD IRF2 

[1] 27
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=8, 14 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 3        3     1.44     1.677       2.4
strata=LOW  5        5     6.56     0.369       2.4

 Chisq= 2.4  on 1 degrees of freedom, p= 0.1 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)      z     p
strataLOW -1.3355    0.2630   0.9231 -1.447 0.148

Likelihood ratio test=2.13  on 1 df, p=0.1445
n= 8, number of events= 8 
   (14 observations deleted due to missingness)


============================

TCGA-COAD PJVK 

[1] 28
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=6, 16 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 2        2     1.62    0.0909     0.141
strata=LOW  4        4     4.38    0.0335     0.141

 Chisq= 0.1  on 1 degrees of freedom, p= 0.7 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)      z     p
strataLOW -0.3503    0.7045   0.9368 -0.374 0.708

Likelihood ratio test=0.14  on 1 df, p=0.712
n= 6, number of events= 6 
   (16 observations deleted due to missingness)


============================

TCGA-COAD CASP5 

[1] 29
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=4, 18 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 3        3     3.75      0.15         3
strata=LOW  1        1     0.25      2.25         3

 Chisq= 3  on 1 degrees of freedom, p= 0.08 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

               coef exp(coef)  se(coef)     z p
strataLOW 2.208e+01 3.903e+09 3.607e+04 0.001 1

Likelihood ratio test=2.77  on 1 df, p=0.09589
n= 4, number of events= 4 
   (18 observations deleted due to missingness)


============================

TCGA-COAD NLRP1 

[1] 30
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=7, 15 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 3        3      3.5    0.0702     0.172
strata=LOW  4        4      3.5    0.0700     0.172

 Chisq= 0.2  on 1 degrees of freedom, p= 0.7 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

            coef exp(coef) se(coef)     z    p
strataLOW 0.3601    1.4335   0.8738 0.412 0.68

Likelihood ratio test=0.18  on 1 df, p=0.6751
n= 7, number of events= 7 
   (15 observations deleted due to missingness)


============================

TCGA-COAD CASP9 

[1] 31
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=9, 13 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 5        5     3.61     0.535      1.02
strata=LOW  4        4     5.39     0.358      1.02

 Chisq= 1  on 1 degrees of freedom, p= 0.3 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)      z     p
strataLOW -0.7416    0.4763   0.7486 -0.991 0.322

Likelihood ratio test=1.02  on 1 df, p=0.3124
n= 9, number of events= 9 
   (13 observations deleted due to missingness)


============================

TCGA-COAD PLCG1 

[1] 32
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=8, 14 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 4        4     4.83     0.142      0.48
strata=LOW  4        4     3.17     0.216      0.48

 Chisq= 0.5  on 1 degrees of freedom, p= 0.5 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

            coef exp(coef) se(coef)     z     p
strataLOW 0.6053    1.8317   0.8851 0.684 0.494

Likelihood ratio test=0.49  on 1 df, p=0.483
n= 8, number of events= 8 
   (14 observations deleted due to missingness)


============================

TCGA-COAD IL18 

[1] 33
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=6, 16 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 3        3     2.43     0.132     0.265
strata=LOW  3        3     3.57     0.090     0.265

 Chisq= 0.3  on 1 degrees of freedom, p= 0.6 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)     z    p
strataLOW -0.4730    0.6231   0.9269 -0.51 0.61

Likelihood ratio test=0.27  on 1 df, p=0.6059
n= 6, number of events= 6 
   (16 observations deleted due to missingness)


============================

TCGA-COAD DPP8 

[1] 34
Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=7, 15 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 4        4      3.8    0.0100    0.0262
strata=LOW  3        3      3.2    0.0119    0.0262

 Chisq= 0  on 1 degrees of freedom, p= 0.9 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

             coef exp(coef) se(coef)      z     p
strataLOW -0.1336    0.8749   0.8268 -0.162 0.872

Likelihood ratio test=0.03  on 1 df, p=0.8716
n= 7, number of events= 7 
   (15 observations deleted due to missingness)


============================

Display the results only for genes where a significant difference in survival has been reported.

significant_genes

[1] "IRF1"  "TIRAP"

num_significant_genes <- length(significant_genes)

if (num_significant_genes > 0) {
  for (i in 1 : num_significant_genes) {
    project <- significant_projects[[i]]
    gene <- significant_genes[[i]]
    
    cat(project, gene, "\n\n")
    gene_df <- construct_gene_df(gene, project)
    
    fit <- compute_surival_fit(gene_df)
    survival <- compute_survival_diff(gene_df)
    cox <- compute_cox(gene_df)
    print(survival)
    cat("\n")
    print(cox)
    print(plot_survival(fit))
    
    cat("\n\n============================\n\n")
  } 
}

TCGA-COAD IRF1 

Warning: Loglik converged before variable  1 ; coefficient may be infinite. 

Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=5, 17 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 3        3     4.35     0.419      4.26
strata=LOW  2        2     0.65     2.804      4.26

 Chisq= 4.3  on 1 degrees of freedom, p= 0.04 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

               coef exp(coef)  se(coef)     z     p
strataLOW 2.232e+01 4.923e+09 3.308e+04 0.001 0.999

Likelihood ratio test=4.61  on 1 df, p=0.03188
n= 5, number of events= 5 
   (17 observations deleted due to missingness)


============================

TCGA-COAD TIRAP 

Call:
survdiff(formula = Surv(overall_survival, deceased) ~ strata, 
    data = gene_df)

n=5, 17 observations deleted due to missingness.

            N Observed Expected (O-E)^2/E (O-E)^2/V
strata=HIGH 2        2     0.65     2.804      4.26
strata=LOW  3        3     4.35     0.419      4.26

 Chisq= 4.3  on 1 degrees of freedom, p= 0.04 

Call:
coxph(formula = Surv(overall_survival, deceased) ~ strata, data = gene_df)

                coef  exp(coef)   se(coef)      z     p
strataLOW -2.232e+01  2.031e-10  3.308e+04 -0.001 0.999

Likelihood ratio test=4.61  on 1 df, p=0.03188
n= 5, number of events= 5 
   (17 observations deleted due to missingness)


============================

---

1. De La Salle University, Manila, Philippines, gonzales.markedward@gmail.com↩︎

2. De La Salle University, Manila, Philippines, anish.shrestha@dlsu.edu.ph↩︎